%%% Example: estimation of a block correlation matrix
clc; close all;

% For selected year there is a T x n matrix of daily returns, X. 
% vB is a K-dimensional vector of block sizes, where K is a number of
% blocks (clusters), and it should be chosen depending on the selected
% block structure (GICS sectors, groups, industries, sub-industries, 
% or the equicorrelation structure). 

% Note: data in matrix X is pre-ordered in such a way that variables
% in the first vB(1) columns correspond to the first cluster,
% variables in the next vB(2) columns are from the second cluster, etc .

% Result (for a pre-selected year and block structure):
% BLKcorr - K x K matrix of estimated block correlations
% BLKcorr_part - K x K matrix of estimated partial block correlations 

% Select year and block structure
year = '2020';           % select a year from 1995 to 2020
block_struct = 'sec';    % can take one of five possible values: 
                         % 'equi' - for equicorrelation structure
                         % 'sec'  - blocks determined by GICS sectors  
                         % 'grp'  - blocks determined by GICS groups 
                         % 'ind'  - blocks determined by GICS industries 
                         % 'sub'  - blocks determined by GICS sub-industries 
                         
if ~any(strcmp(sprintfc('%d',(1995:2020)'),year))
    fprintf ('Error : wrong year \n'); return
end
if ~any(strcmp({'equi','sec','grp','ind','sub'},block_struct))
    fprintf ('Error : wrong block structure name \n'); return
end
 
% Load X and vB
X = dlmread(['crisp_returns_',year,'.csv'],',',1,0); 
X = X(:,2:end); [T,n] = size(X);
if strcmp(block_struct,'equi')
    vB = n;
else
    vB = dlmread(['gics_',block_struct,'_',year,'.csv'],',',1,1);
end

% Estimation of the canonical form
X = X./sqrt(mean(X.^2));   % Standardized returns
[A_hat,Lambda_hat,logLik,BIC] = MLE_block_cov(X,vB);

% Block correlation estimates
K = length(vB); BLKcorr = zeros(K); 
for i = 1:K
    for j = 1:i
        if j == i   % within block correlations (only if vB(i) > 1)
            if vB(i) > 1
                BLKcorr(i,i) = (A_hat(i,i)-1)/(vB(i)-1);   % or (A_hat(i,i)-Lambda_hat(i))/vB(i)
            else
                BLKcorr(i,i) = NaN; 
            end
        else        % between blocks correlations 
            BLKcorr(i,j) = A_hat(i,j)/sqrt(vB(i)*vB(j));
            BLKcorr(j,i) = BLKcorr(i,j);
        end
    end
end

% Matrix of conditional correlations
BLKcorr_part = zeros(K); 
for i = 1:K
    for j = 1:i
        if j == i   % cond. correlation for a pair of stocks from the same block
            A_cond = A_hat(i,i) - ...
                A_hat(i,setdiff(1:K,i))*(A_hat(setdiff(1:K,i),setdiff(1:K,i))\A_hat(setdiff(1:K,i),i));
            if vB(i) > 1
                BLKcorr_part(i,i) = (A_cond(1,1)-Lambda_hat(i))/(A_cond(1,1)+(vB(i)-1)*Lambda_hat(i));
            else
                BLKcorr_part(i,i) = NaN;
            end
        else        % cond. correlation for a pair of stocks from two distinct blocks
            A_cond = A_hat([i,j],[i,j]) - ...
                A_hat([i,j],setdiff(1:K,[i,j]))*(A_hat(setdiff(1:K,[i,j]),setdiff(1:K,[i,j]))\A_hat(setdiff(1:K,[i,j]),[i,j]));
            if vB(i) > 1
                di = (vB(i)-1)*Lambda_hat(i)+A_cond(1,1);
            else
                di = A_cond(1,1);
            end
            if vB(j) > 1
                dj = (vB(j)-1)*Lambda_hat(j)+A_cond(2,2);
            else
                dj = A_cond(2,2);
            end
            BLKcorr_part(i,j) = A_cond(1,2)/sqrt(di*dj);
            BLKcorr_part(j,i) = BLKcorr_part(i,j);
        end
    end
end
